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Abstract 

We use a series of cosmological A^-body simulations and various analytic models 
to study the evolution of the matter power spectrum in real space in a A Cold Dark 
Matter universe. We compare the results of A^-body simulations against three an- 
alytical model predictions; standard perturbation theory, renormalized perturbation 
theory, and the closure approximation. We include the effects from finite simulation 
box size in the comparison. We determine the values of the maximum wavenumbers, 
k^l^^ and /cgo^, below which the analytic models and the simulation results agree to 
within 1 and 3 percent, respectively. We then provide a simple empirical function 
which describes the convergence regime determined by comparison between our sim- 
ulations and the analytical models. We find that if we use the Fourier modes within 
the convergence regime alone, the characteristic scale of baryon acoustic oscillations 
can be determined within 1% accuracy from future surveys with a volume of a few 
h~^Gpc^ at 2 ~ 1 or 2 ~ 3 in the absence of any systematic distortion of the power 
spectrum. 

Key words: cosmology: large-scale structure of universe — methods: statistical 
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1. Introduction 



The nature of dark energy remains one of the most fundamental questions in physics 
and cosmology. It is accessible only through astronomical observations, and a number of large 
galaxy redshift surveys and weak-lensing observations are expected to yield tight constraints 
on the dark energy equation of state parameter, w^e =Pde/pde- The use of baryon acoustic 
oscillations (hereafter BAOs) is a promising tool to determine wde- The characteristic scale 
of BAOs can be used as a robust standard ruler, which helps us to reconstruct the expansion 
history of the universe. 

Substantial improvements in theories as well as observations are required to constrain 
wbe within an accuracy of a few percent, which is the goal of next generation BAO surveys. 
Nishimichi et al. (2007) showed that the measurements of the angular diameter distance and 
of the Hubble parameter must be very accurate, with errors less than 1%, in order to constrain 
Wbe at a 5% level at 2; ~ 1 in the special case that the other cosmological parameters are fixed. 
Clearly, theoretical models are required to predict the matter power spectrum with even better 
accuracy. 

Recently, significant progress has been made in developing analytic models of the matter 
power spectrum by extending conventional perturbation theory (e.g., Crocce & Scoccimarro 
2006a, b, 2008; Matarresse & Pietroni 2007, 2008; Valageas 2007; Taruya & Hiramatsu 2008; 
Matsubara 2008). Cosmological A^-body simulations are often used as a reference to calibrate 
these models. It is not clear, however, whether or not A^-body simulations provide sufficiently 
accurate results at the required percent level even at weakly nonlinear scales where perturbation 
theory is generally supposed to be reasonably accurate. 

Indeed a variety of systematic effects needs to be considered in interpreting the results of 
A^-body simulations. It is well known that discreteness of a density field in A^-body simulations 
causes spurious behaviour of the power spectrum at length scales comparable to or smaller than 
the mean inter-particle separation (Melott et al. 1997; Splinter et al. 1998; Hamana, Yoshida & 
Suto 2002; Baertschiger, Joyce & Sylos Labini 2002; Marcos et al. 2006; Joyce & Marcos 2007a, 
b; Joyce, Marcos & Baertschiger 2008; Romeo et al. 2008). Adopting a finite-size simulation box 
and imposing a periodic boundary condition also systematically biases the growth of density 
fluctuations at almost all the length scales (Seto 1999). Recently, Takahashi et al. (2008) 
examined the effect of the finite number of Fourier modes in a single A^-body realization due 
to a finite box size. They found that the evolution of the power spectrum measured from each 
realization deviates from the prediction of perturbation theory by more than a few percent 
even at very large scales (e.g. k ^ O.l/iMpc"^) and that this deviation results from having 
only a finite number of modes, which can be explained by including second-order perturbation 
contributions. This implies that a sub-percent level accuracy both in theories and in simulations 
is very demanding and we need to explore this in detail. 
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In the present paper, we consider a specialized case of the matter power spectrum in 
real space. We critically compare the evolution of the matter density fluctuations in the weakly 
nonlinear regime in cosmological simulations to that of theoretical predictions. We measure the 
matter power spectra using the simulation outputs and compare them with analytic models. In 
doing so we correct for the finite-mode effect pointed out by Takahashi et al. (2008). Unlike the 
conventional approach, we do not regard the results of A^-body simulations as perfect calibrators 
of analytical models. We determine ranges of wavenumbers where both theories and A^-body 
simulations agree within a given accuracy. 

The rest of the present paper is organized as follows. Section 2 briefly outlines the 
analytic models that we examine for the nonlinear matter power spectrum. Section 3 describes 
our simulation details. Our methods to measure the matter power spectrum and to correct for 
the finite volume effect are shown in section 4. We show our results of A^-body simulations 
and determine the convergence regime in wavenumber in section 5. We also discuss the phase 
information of BAOs and predict the parameter forecast from the modes in this convergence 
regime alone in section 6. Finally section 7 summarizes the present paper. The convergence 
tests of our A^-body simulations are discussed in the appendix. 

2. Analytic Nonlinear Models 

Various analytic models have been proposed to account for nonlinear evolution of the 
matter power spectra in real space. In this paper, we will compare A^-body simulations with 
three different analytic models, specifically focusing on the leading-order contributions. The 
models include standard perturbation theory (SPT; e.g., Bernardeau et al. 2002), renormalized 
perturbation theory (RPT; Crocce & Scoccimarro 2006a, b, 2008), and the closure approxima- 
tion (CLA; Taruya & Hiramatsu 2008). 

SPT is a straightforward expansion of the fluid equations around their linear solution, 
assuming that fluctuation amplitudes are small. Schematically, the expansion is written 

P{k, z) = P'^ik, z) + Pi-^°°P(A:, ^) + ■ • • , (1) 

where P^{k,z) is the linear power spectrum, which grows as oc D'^{z) with linear growth rate, 
D_^.{z). The second term, P^~^°°'^{k, z), is the one-loop correction to the power spectrum, 
which represents the contributions coming from the second-order and third-order solutions. 
The one-loop correction is roughly proportional to P^Aq with Aq = k^P^{k,z)/ (2tt^), and this 
term subsequently exceeds the linear term at lower redshift and smaller scales. The explicit 
expressions for the one-loop correction P^~^°°'^{k,z) may be found in the literature (e.g., Makino, 
Sasaki & Suto 1992; Jain & Bertschinger 1994; Scoccimarro & Frieman 1996; Jeong & Komatsu 
2006; Nishimichi et al. 2007). 

Both RPT and CLA are constructed from the renormalized expression of a perturba- 
tion series obtained from SPT, by introducing non-perturbative statistical quantities such as 
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the propagator and the vertex function. The resultant renormalized expressions for the power 
spectrum are given as an infinite series of irreducible diagrams composed of the propagator, the 
vertex function, and the power spectra themselves in a fully non-perturbative way. Then, em- 
ploying the Born approximation, RPT of Crocce & Scoccimarro (2008) perturbatively evaluates 
the renormalized diagrams under the tree- level approximation of vertex functions. In contrast, 
CLA proposed by Taruya & Hiramatsu (2008) first truncates the renormalized expansion at 
the one-loop order under the tree-level approximation of vertex functions, and obtains a closed 
set of equations for the power spectrum and the propagator. In practice, they further apply 
the Born approximation to the truncated diagram. As a result, the leading-order calculations 
of the power spectra in CLA reduce to the same analytic expression as in the case of RPT, 
leading to 

P(fc, z) = G\k, z)P^{k, z,^0 + P^c'°°P(fc, z) + ---. (2) 

Here, the function G{k,z) is the propagator, and the term Pl^°°^{k,z) represents the corrections 
generated by mode-mode coupling at smaller scales, constructed from the one-loop diagram. In 
the above expression, the only difference between CLA and RPT is the asymptotic behavior of 
the propagator. The explicit expressions for the propagator together with the mode-coupling 
term are summarized in Crocce & Scoccimarro (2006b, 2008) and Taruya & Hiramatsu (2008). 

We note here that the most important difference between SPT and RPT/CLA is the 
convergence properties dictated by the propagator. Similar to P^~^°°^{k,z) of SPT, the term 
Pl^°°'^{k,z) is roughly proportional to P^Aq, but it additionally contains propagators, leading 
to a large suppression at high k. Therefore, to the one-loop level, RPT and CLA are expected 
to be most accurate at low k (below the propagator damping scale), but the predictions at low 
k are expected to be much accurate for RPT/CLA compared to SPT. 

3. Simulations 

In this section we describe our fiducial simulations. We also run some simulations with 
different settings (initial conditions, integrators and box sizes), which are discussed in the 
appendix. 

We begin with an explanation of our initial conditions for the A^-body simulations. We 
compute the linear power spectrum at an initial redshift Zini = 31 using CAMB (Lewis, Challinor 
& Lasenby 2000). For all the simulations in the present paper, we adopt the standard ACDM 
model with the cosmological parameters from the WMAP3 results (Spergel et al. 2007), Qm = 
0.234, = 0.766, Qb = 0.0414, h = 0.734, dg = 0.76 and n, = 0.961, which are the current 
matter, cosmological constant, baryon densities in units of the critical density, the Hubble 
constant in units of 100km s~^Mpc~^, density fluctuation amplitude smoothed with a top-hat 
filter of radius 8h~^Mpc, and the scalar spectral index, respectively (see table 1). We then 
generate a linear overdensity field in Fourier space assuming Gaussianity; the amplitude follows 
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Table 1. Adopted cosmological and simulation parameters of Gadget-2 for our fiducial runs. 



cosmological 








h 




ns 


value 


0.234 


0.766 


0.175 


0.734 


0.76 


0.961 


simulation 


box size 


# of particles 




# of PM grids 


softening length 




value 


lOOO/i-^Mpc 


5123 


31 


10243 


O.l/i-^Mpc 


4 



Gaussian sampling particle 
n-th realization realization 



N-body 
simulation 



square and 
average over 
k and n 



5b (z, ,) — ^ (5N-'""*y(z. .) 

k.nV ini' k,n \ ini' 



perturbation theory 



5N-body(z) 



pN-body(k^,z) 



square and 
average over 
k and n 



linear theory 



^ PKkpz) 

Fig. 1. A flow chart to illustrate our methodology to correct for the effect of finite box size. 

the Rayleigh distribution, and the phase is uniformly distributed. We employ 512^ dark matter 
particles within a cube of 1000/i~^Mpc in each side. We displace these particles from regular 
grid pre-initial positions using the second-order Lagrangian perturbation theory (2LPT; e.g., 
Crocce, Pueblas & Scoccimarro 2006). 

We next describe the time integration scheme. We adopt a publicly available parallel 
cosmological A^-body solver Gadget2 (Springel 2005). The number of meshes used in the 
particle-mesh computation is 1024^. We adopt a softening length of 0.1/i~^Mpc for the tree 
forces. We select three output redshifts; 2 = 3,1 and 0, where we measure the power spectrum 
(see section 4 for more details). 

In the appendix, we discuss the convergence of our simulations against different initial 
conditions, box sizes, and codes. We found that the current setup provides a convergent result 
within 1% of the amplitude of the power spectrum at our scales of interest {k ^ 0.4/iMpc~^). 
We use this setup as our fiducial model, and run 4 different realizations for this model. 
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Table 2. Notations for various density contrasts and power spectra. 





meaning 


description 




input linear power spectrum 


text m gz 




density contrast Gaussian-sampled from P^{k,z) 
for the n-th realization 


text in §3 


s-W-body / \ 


density contrast realized by particles using 2LPT displacement 
and evolved by A'^-body simulation 


text in §3 


pJV-body^L. _\ 


pOWci opeCbiUlll Ui Lllc t-Lli WdVcliUlllUei Ulll toLillldLtt-l ii Olll iV-UOLiy 

simulations taking average over finite modes and realizations 






^kni^) evolved by perturbation theory 


Eq.(ll) 




same as P^-^°'^y{ki,z) but calculated from Sll{z), not S^~^°'^y{z) 


Eq.(lO) 


pAT-body/, N 
^corrected V^-ii^y 


pTV-body^^^ corrcctcd for the effect of finite volume 


Eq.(12) 



4. Power Spectrum Analysis 

4.1. Matter density field and power spectrum in N-body simulations 

Here we briefly describe the notation for various density fields and power spectra which 
we use throughout the paper (see figure 1 and table 2). 

We denote the Gaussian-sampled density field from P^{k,Zi^i) by 5^,„(2ini), where the 



Af-body / 
k,n \^ 

is different from 6, 



the density field 

(^ini) ) 



Af-body 
k,n 



subscript n stands for the n-th realization. We then denote by 6 
realizedhj particles using 2LPT displacement. Note that S^^{. 
the former is strictly Gaussian, whereas the latter slightly deviates from Gaussian because 
of the nonlinearity and discreteness in the process of displacement. The measured density 



contrast from the simulation output at redshift z is also expressed as ^j^. 

5, 



TV-body 



z). In measuring 



A^-body / 
k,n \ 



from the simulation output, we first assign particles onto a 1024"^ mesh using Cloud- 
in-Cells interpolation (Hockney & Eastwood 1981). We then use a FFT to calculate density 
contrasts in Fourier space, and divide each Fourier mode by the Fourier transform of the window 
function (Angulo et al. 2008). We made sure that the details of the interpolation scheme and 
the number of mesh points do not affect the results significantly in the scales of our interest 
{k ^ OAh Mpc ~^). See also Jing (2005) for a discussion of the effects of aliasing. 



We square 6^^ ^"'^^{z 

pAf-body^^.^^)^ 

Jv-i =^ 



and take an average over realizations and modes: 

1 



A^-body / 



k,?i 
1 



AT" 



jymode jyr 



<|k|<fc™ n=l 



E 



(3) 
(4) 



where A"^"'^^ and A'^'^" are the numbers of modes in the i-th wavenumber bin and the number of 
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realizations, and /cf ™ and kf^^^ are the minimum and the maximum wavenumber, respectively. 
Note that we use to denote the average over modes in the i-th wavenumber bin and over 
realizations: this average is not equivalent to the true ensemble average, and the difference 
corresponds to the finiteness of the simulated volume (or number of modes). In what follows, 
we adopt equally-spaced bins with width Ak = 0.005/iMpc~^. 

Finally, the standard errors of the averaged power spectra of equation (3) can be esti- 
mated by 



error of P^-^°'^y(A;i 



(5) 



jymodc jyrun 

Note that this value indicates the uncertainty in the estimation of the central value, not the 
variance among the modes in each bin. 

4-2. Corrections to the power spectrum 

The matter power spectrum measured from simulations deviates from the prediction for 
the ideal ensemble average, which can be obtained only in the limit of an infinite number of 
realizations or an infinite box size. This deviation is actually important in interpreting the 
results of A^-body simulations as shown by Takahashi et al. (2008); the matter power spectrum 
measured from their A^-body simulations does not agree with the predictions of linear theory 
nor SPT even at very large scales (e.g., k ^ 0.1/iMpc~^). They examined the effect of a finite 
box size (hence a finite number of modes) and showed that the finite-mode effect is actually 
responsible for the anomalous growth rate. Here we briefly summarize their formulation of the 
correction. 

We follow the standard notation used in cosmological perturbation theory [see 
Bernardeau et al. (2002) for a review]. Let us expand the density perturbations in /c-space 
for the n-th A^-body realization as 

Cn-'°'^W=5tW+CW + -- (6) 

Here the second-order term is expressed as a sum of contributions from mode couplings between 
two modes: 

C(^) = E^^'nP,k-p;.)5;;,„(z)5Lp,nW, (7) 
p 

where the symmetrized second-order kernel F^'^^ is expressed as 

1 r 1 x ■ y 

e{z) ^ (9) 

In practice, the time dependence of the kernel function is very weak, and thus we simply set 
e = 3/7 in what follows. The power spectrum up to the third order in S]^^{z) averaged over 
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F(^)(x,y; .) = i [1 + eiz)] + f - + ^ ) + ^ [1 - e{z)] (8) 



modes in the i-th wavenumber bin and over realizations is 

2\ 



23? 



(10) 

(11) 



where 9?[...] stands for the real part of a complex number. Though the second term should 
vanish for an ensemble average over infinite modes, it does not vanish exactly if the number 
of Fourier modes is finite. The first term grows as oc D'^{z) while the second term grows as 
oc D'^{z), thus the second term becomes important at late time (i.e., at low redshifts). 

The method to correct for the deviation from the ideal ensemble average is to multiply 
P^-^°'^y{ki,z) by the ratio of P^^ih^z) and P^ih^z): 

Pco;re:ttih,z) ^ P''-^°''y{h,z) X P^ (h, z) / P^^ {h, z) . (12) 

The individual random nature of each A^-body run in both P^~^°'^y{ki, z) and P^^{ki,z) is 
weakened by this procedure [equation (12)] as long as the predictions of perturbation theory 
are sufficiently accurate. This method does not bias the corrected value of the power spectrum 
even if the perturbation theory breaks down, since P^"^{ki,z) approaches P^{ki,z) in the limit 
of infinite number of Fourier modes. 

As in equation (5), the standard errors of equation (12) can be estimated as 



C nA'^— body 

error of -r. 



corrected 



(,ki,z) 



■Af— body / 
k,n [ 



T 2\ 



The average of 



A'— body / 
k,n \ 



is P^'^°'^y{ki,z) while that of 



n jymodc 



(13) 



is P {ki,z), and thus we 



multiply the ratio in the second term in the numerator to adjust the mean values. 
5. Results 

5.1. Comparison between N-body simulations and analytic models 

As mentioned before, the accuracy of A^-body simulations themselves is not perfect and 
we do not regard them as perfect calibrators of theoretical models. Instead we compare the 
power spectrum of our simulations with those from theoretical predictions aiming at deter- 
mining the reliable range of wavenumbers in which both simulations and theoretical models 
agree. 

We first compare the averaged power spectrum over the four realizations without any 
corrections to the theoretical predictions. The left panel of figure 2 plots the fractional differ- 
ences from the linear power spectrum, P'^™(A;,z), computed from no- wiggle formula of Eisenstein 
& Hu (1998). The errorbars indicate the standard errors of the estimated mean value [equation 
(5)]. To clarify the differences between the A^-body results and theoretical predictions, we also 
plot the residuals from RPT in the right panel of figure 2. 
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Fig. 2. Comparison of simulation power spectra and analytical predictions: Left: Power spec- 
tra of our simulations before correction, normalized by the no- wiggles formula; top 2 = 3, mid- 
dle z = 1, bottom z ~ 0. The errorbars show the standard errors [equation (5)]. Lines are the- 
oretical predictions described in section 2; dotted: standard perturbation theory (SPT), dot-dashed: 
renormalized perturbation theory (RPT), dashed: closure approximation (CLA), solid: linear the- 
ory (LIN). Right: Same as the left panel but we plot the difference from the RPT prediction. 

Since large errorbars aX k ^ 0.1/iMpc"^ are expected to come mostly from the finite 
volume effect, we next correct for this. In figure 3, we plot the power spectra, based on the 
procedure in § 4.2, but we truncate the expansion of equation (11) at the first term {left: the 
fractional difference from the no- wiggles formula, right residuals from RPT prediction). The 
errorbars become significantly smaller compared with those in figure 2, because the finiteness 
effect is reduced by our methodology. Nevertheless, the results of A^-body simulation still 
exhibit a few percent error, even after the correction. 

Next we include the second term of equation (11) for the correction, which is coming 
from the mode couplings among finite modes. Figure 4 shows the results for both the fractional 
difference from the no-wiggles formula (left), and the residuals from RPT predictions (right). 
Now the size of errorbars are further reduced to the sub-percent level. 

All the theoretical predictions plotted in figure 4 and A^-body simulations agree with 
each other well within the errorbars at large scales up to some wavenumbers (we will determine 
the range of convergence in the next subsection). Among the four theoretical predictions, linear 
theory deviates from the rest at the smallest wavenumber. The range of agreement is wider for 
SPT than linear theory, as SPT includes the leading-order contribution of nonlinear growth. 
RPT and CLA seem to show better agreement with A^-body simulations compared with SPT 
although all of the three nonlinear models include their own leading-order nonlinear corrections. 
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Fig. 3. Same as figure 2, but we correct for the effect of finite volume. We trun- 
cate the expansion of equation (11) at the first term. The errorbars show equation (13). 
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Fig. 4. Same as figure 2, but we correct for the effect of finite vohime includ- 
ing the second term of equation (11). We also show the 1% limit wavenum- 
bers, fcj*^, for LIN, SPT and RPT/CLA by vertical arrows (from left to right). 

This difference in the agreement ranges corresponds to their different convergence properties; 
RPT and CLA show better convergence at scales where nonlinearity is very weak. 
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5.2. Convergence regime in wavenumber 

We are now able to quantitatively estimate the convergence regime of wavenumber where 
theories and A^-body simulations agree. We define two characteristic wavenumbers, k^^^ and 
/cg™ , such that the results of A^-body simulations and theoretical prediction agree within 1% at 
k < k^^^ and within 3% at k < /cg™, respectively. We confirmed that if we add 1% (3%) Gaussian 
errors on the power spectra binned by Ak = 0.005/iMpc~^ that leads to a ~ 2% (~ 6%) error 
on the dark energy equation of state parameter wde at z = 1 using the phase information of 
BAOs, which is roughly the goal of upcoming (ongoing) surveys. 

Before determining the wavenumbers, we briefiy summarize three convergence criteria 
previously introduced in the literature. The first one is 

k^P^_^{k,z) 

introduced by Jeong & Komatsu (2006) as a 1% convergence regime of SPT. The second one 
is introduced by Sefusatti & Komatsu (2007): 

a^Rmin, z)= [ d'kW'iqRmin)P''{q, z) = 0.25, (15) 



A^(fe,^)^ o r ^ <0A (14) 



TT 

^max = -Z-^ 5 (16) 

where H^(q'-Rmin) denotes the Fourier transform of a top-hat window function with radius -Rmin- 
The third one is proposed by Matsubara (2008): 

A;V2 = — P'^{q,z)dq<Q.2h, (17) 

where cr^ is the one-dimensional linear velocity dispersion. This gives a scale where nonlinearity 
becomes important. Note that Crocce & Scoccimarro (2006) presented similar criteria. 

We show these criteria by dotted lines against redshift in figure 5. The dotted curve 
labeled with "JK06" represents the k\% determined from equation (14), and the /cmax of equation 
(16) is shown by "SK07", and the nonlinear scale of equation (17) is shown by "M08". Among 
these three criteria, the first two have similar redshift dependence. This is because they are 
based on a local value of the power spectrum or its Fourier transform convolved with a window 
function. The other criterion of Matsubara (2008) refers to the integrated value of the power 
spectrum over scale and thus this shows a rather mild redshift dependence of the maximum 
wavenumber. 

The values of k\% and k\% can be directly read off from figure 4, which are shown in table 
3 and also in figure 5 {k\% by filled symbols and fcg™ by open symbols, squares: RPT/CLA, 
triangles: SPT, circles: LIN). We empirically find that a modified version of equation (17): 

^^l'p\q,z)dq<C, (18) 

well reproduces the results in a unified fashion. The constant C in the right-hand-side depends 
on the choice of theoretical model and the threshold value (1% or 3% in this paper). We find 
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Table 3. The values of fcj'™ and fc^'™ read fr om figure 4 and the corresponding eonstant value in equation (18). 





z = 3 


z = l 


z = 


C-1% 


2 = 3 


z = l 


z = 




RPT/CLA 


0.3 


0.18 


0.12 


0.35 


0.36 


0.20 


0.14 


0.5 


SPT 


0.22 


0.13 


0.08 


0.18 


0.29 


0.16 


0.11 


0.3 


LIN 


0.13 


0.09 


0.06 


0.06 


0.19 


0.12 


0.08 


0.13 



that = 0.35, C§^^ = 0.18 and C^™ = 0.06 well reproduce the 1% agreement limit of 

RPT/CLA, SPT and linear theory from figure 4, respectively. We plot equation (18) with these 
constant values as solid lines in figure 5 and also as vertical arrows in figure 4. Similarly we 
find that the corresponding 3% limits are given by Cfc^'^'^^^ = 0.5, C§%^ = 0.3 and = 0.13 
for the three theoretical predictions (dashed lines in figure 5). 

The convergence regimes of our criteria [equation (18)] are narrower than those previ- 
ously proposed. Note, however, that the criteria (18) with the above constant values reasonably 
describe the valid range of analytic models even at higher redshifts {z = 7 and 15). One of the 
possible sources of the difference between our criterion and the previous ones is that they use 
simulations with smaller box sizes to achieve better resolution, which enhances a systematic 
error due to the finiteness of box size. More details are discussed in the appendix. 

6. Implications for the Phases of BAOs 

6.1. Extraction of the phases of BAOs from the nonlinear power spectra 

So far we have focused on the amplitudes of the matter power spectra, but the phase 
information of BAOs is also important. This is imprinted mainly in the leading term which 
is the product of the linear power spectrum and the propagator in RPT and CLA [see equa- 
tion (2)]. Since the other contributions, the mode-coupling terms, are fairly smooth functions 
of wavenumber, the phase information of BAOs is almost erased by the convolution of the 
smoothed kernel and the propagator ^. Therefore almost all of the BAO information comes 
from the zeroth-order term in RPT/CLA. We thus expect that the phase of BAOs predicted 
by A^-body simulations and theories agree with each other for a wider range of wavenumbers 
compared to the amplitudes since the amplitudes are sensitive to higher-order corrections. 

We test the convergence of the predicted phases of BAOs as follows. We first construct a 
smooth reference power spectrum using the spline technique developed by Percival et al. (2007) 
and Nishimichi et al. (2007). We adopt a basis spline (B-spline) fitting function as the reference 
spectrum, P^-^p"'^^(A;), with the break points at k}^''''-'^ = O.OOl/iMpc-^ and kf""^^ = (0.05j - 
0.025)/iMpc"i where j is a positive integer. The data points to be fitted are set at ki of 

^ The mode-coupling term has very small wiggly feature which comes from BAOs. This leads to a sub-percent 
level shift of the positions of "nodes" . See Crocce & Scoccimarro (2008) for more details. 
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Fig. 5. Upper limit of reliable wavenmiibers, fc'™ and fcg™, described in the text. Symbols show the 
values read from figure 4. circles: linear theory, triangles: SPT, squares: RPT/CLA. Filled sym- 
bols correspond to fcj'o™, while open ones represent ^3™. The three solid lines plot equation (18): 
fci'^ for linear theory (C = 0.06), SPT (C = 0.18) and RPT/CLA (C = 0.35) from left to right, 
and the dashed lines are corresponding fcg™ (C = 0.13, 0.3 and 0.5, respectively). We also show 
the nonlinear wavenumbers proposed by Jeong & Komatsu (2006), Sefusatti & Komatsu (2007) and 
Matsubara (2008) as dotted lines with JK06, SK07 and M08. The two shaded regions show the rcdshift 
range planed by WFMOS survey (with minimum wavenumbers, li^ jV^I'^ . Here V is the survey volume.). 

equation (4) with a bin width of iSk = O.OOS/iMpc"^ for both the A^-body simulations results 
and the theoretical predictions. We use the A^-body power spectrum, -PcOTrected(^' which 
the effect of finite volume is corrected using perturbation theory. We assign equal weights to 
all data points. A smooth reference spectrum is computed for the average of the simulations 
and the theoretical models. 

We then divide the model and A^-body power spectra by their respective smooth refer- 
ence power spectra, P^~'^p'"^'^(/c), described above, which are shown in figure 6. The symbols 
and lines have the same meanings as in figure 4. Among the three nonlinear models, RPT 
and CLA lie almost on top of each other, and they are indistinguishable. These two models 
show better agreement with A^-body simulations than SPT, whose predicted strong damping of 
oscillations incorrectly leads to a phase reversal {k ^ 0.25/iMpc~^ at 2; = 1, A; ^ 0.16/iMpc~^ at 
z = 0). Compared with the vertical arrows which represent the values of /c[™ determined by the 
amplitudes of the power spectra, the convergence regimes of BAO phases predicted by analytic 
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models and iV-body simulations are quite a bit wider. We thus conclude that the phases of 
BAOs are potentially more robust and more useful in the actual analysis of BAO surveys. 
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Fig. 6. Power spectra divided by tlieir B-spline fits described in tlie text. 
Symbols, lines and vertical arrows are analogous of those in figure 4. 

6.2. Recovery of the BAO scale from WFMOS survey 

With our criteria for trustable ranges of simulations and analytic models, we are able to 
discuss future cosmological constraints using BAOs. A number of simulation papers attempted 
to present parameter forecasting of the dark energy equation of state parameter, wbe =Pde/pde, 
using BAOs as a standard ruler (i.e., Meiksin, White & Peacock 1999; Springel et al. 2005; 
Angulo et al. 2005, 2008; Eisenstein, Seo & White 2007; Eisenstein et al. 2007; Seo & Eisenstein 
2003, 2005; Seo et al. 2008; Smith, Scoccimarro & Sheth 2007, 2008; White 2005; Huff et al. 
2007; Jeong & Komatsu 2006, 2008). We follow a similar procedure, taking into account the 
reliable range of wavenumbers in the prediction of the BAO scale. 

We consider the Wide-Field Multi-Object Spectrograph (c.f., Bassett, Nichol & 
Eisenstein 2005, WFMOS) as a specific example in the present analysis. We construct the 
template power ratio, P/P^~^p^™°(/c), using RPT and linear theory as described above. We cre- 
ate mock power spectra by adding Gaussian random errors (Feldman, Kaiser & Peacock 1994): 
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(19) 



' P{k) \ _ APjk) _ I 2 

We neglect the errors arising from the construction process of P^~'^p^™°(A;). In the above ex- 
pression, iV™°'^'^^ is the number of modes in the wavenumber bin k ^ k + Ak: 

^modes ^ (20) 

where V denotes the survey volume, and Ug in equation (19) is the mean number density of 
galaxies. We set Ak to be very small (0.001/iMpc~^) so as not to affect the results and use the 
nonlinear matter power spectrum of Smith et al. (2003) for the calculation of the shot noise 
contribution, [ngP{k)]~^ , just for simphcity. 

We fit these mock data to the corresponding theoretical template P/P^~'^P^™*^(a/i;) with 
a single free parameter a, the ratio of the true distance scale to the assumed one. The response 
of the observed power spectrum, Pohs{k,z), to the parameter a can be written as (e.g., Jeong 
& Komatsu 2008) 

Pohs{k,z) =a~^Ptrnc{ak,z); a = ^^^^, (21) 

-t^V,fid 

where Ptrue(^5^) is the true power spectrum and Dy^z) oc [D\(z)/H(z)Y^^, where D^iz) is 
the angular diameter distance (Eisenstein et al. 2005). -Dv,truc(-2) represents the true distance, 
while Dyfid{z) is calculated assuming the fiducial cosmology. In fitting the templates, we 
fix the value of the minimum wavenumber, fcmin; to be 0.02/iMpc~^, while the value of the 
maximum wavenumber, fcmax, is varied. The assumed parameters for the z ~ 1 (z ~ 3) WFMOS 
survey are taken from Glazebrook et al. (2005): V = 4h~^Gpc^ {lh~^Gpc^) for the volume, 
Ug = 5x 10~^/i^Mpc~^ (5 X 10~^/i^Mpc~^) for the galaxy number density, and b = 2.0 (3.5) for 
the linear bias parameter. 

We employ a Markov chain Monte Carlo fitting routine to determine the uncertainty in 
the scale shift, cTq,, adopting RPT for the fiducial spectrum. We plan to do a similar fitting with 
more free parameters in future work. Figure 7 show the uncertainty in the scale shift, (7 , as a 
function of /cmax (solid lines). We also plot the results when we do not include shot noise and 
nonlinear degradation of BAOs (linear template) by the dashed lines. As a reference, vertical 
solid and dashed arrows indicate k^^!^ for RPT and linear theory, respectively. The difference 
of solid and dashed lines correspond to the impact of nonlinearity and shot noise, i.e., the 
information of phase of BAOs is lost due to these effects. The resulting error of a is 0.7(1.0)% 
at z = 1(3) when we adopt k^^^ of RPT/CLA as /Cmax- These values are slightly improved and 
become 0.6(1.0)% if we increase fcmax to 0.4;iMpc-^ We can put a strong constraint on a even 
if we do not include the modes at A; > k^^^ of RPT in the analysis. In reality, however, we need 
to take account of galaxy clustering in redshift space in addition to what we did in the current 
paper, which will be presented in future work. 
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Fig. 7. Estimated constraints on the scale shift parameter, a, for both the wide (z ^ 1, red) and the 
deep (z ~ 3, blue) components of the proposed WFMOS survey. The sohd lines show the results 
when we use RPT as the template spectrum, while the dashed ones correspond to an ideal case of 
no shot noise and no nonlinear degradation of BAOs (using the linear theory template). We fix the 
minimum wavenumber to be used in the fitting, k^in = 0.02/iMpc~^, and the results are plotted as 
a function of the maximum wavenumber. The assumed survey parameters are written in the figure. 

7. Summary 

We have carried out a systematic comparison of power spectra calculated from A^-body 
simulations and various theoretical models. We correct for the effect of finite volume of A^- 
body simulations using the perturbation-theory-based methodology developed in Takahashi 
et al. (2008). 

We find that our simulations agree with all the theoretical models used in the present 
paper at large length scales {k ^ O.l/iMpc"^). At smaller length scales where nonlinearity 
is mild, our simulations, RPT (Crocce & Scoccimarro 2006a, b, 2008) and CLA (Taruya & 
Hiramatsu 2008) agree with 1% and 3% accuracies up to k^^% and fcg™ as presented in figure 5. 
These convergence regimes depend on the redshift and can be explained by a simple empirical 
formula of equation (18). We also showed that the phase information of BAOs extracted from 
power spectra using B-spline fitting is more robust than the amplitudes: predicted phases of 
nonlinear theories and A^-body simulations agree well in a wider range of wavenumbers than 
the amplitudes of the power spectrum. 
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The currently achieved accuracies of both theoretical predictions and simulations are 
sufficient to interpret the data from future surveys like WFMOS: one can put a tight constraint 
on the BAO scale (cTq, ^1%) using modes within our convergence regimes. There is the potential 
to extract more than simply the BAO scale from the power spectrum, such as the sum of the 
neutrino masses, if one is able to exploit smaller scale information (e.g., Saito, Takada & 
Taruya 2008). To this end, it is important to consider higher-order corrections (i.e., 2-loop 
terms) on the theory end and higher resolution simulations. It is also important to investigate 
the accuracy of the velocity field in A^-body simulations to model accurately in redshift space. 
These are future investigations and a study in redshift space is now in progress. 

We also checked the convergence of the matter power spectrum in A^-body simulations 
changing the initial conditions, A^-body solvers, and box sizes. Among these we showed that A^- 
body simulations with small box sizes (^ 500/i~^Mpc) suffer from systematic enhancements of 
the matter power spectra at BAO scales (e.g., Springel et al. 2005; Seo & Eisenstein 2005; Jeong 
& Komatsu 2006). 
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Appendix 1. Convergence check of different A/^-body codes and dependence on 
simulation parameters 

A. 1.1. Initial Condition Generator 

We test two methods for generating cosmological initial conditions. One is the Zel'dovich 
approximation (Zel'dovich 1970; hereafter ZA) commonly used for cosmological simulations, 
and the other is second-order Lagrangian perturbation theory (2LPT; e.g., Crocce, Pueblas & 
Scoccimarro 2006). Starting from the same linear density field, we generate two initial 
conditions at z^^i = 31 from A-body runs using ZA and 2LPT. Then we compare the power 
spectra evolved with Gadget2. 

Figure 8 shows the ratio of the two matter power spectra from simulations at various 
output redshifts. The overall trend that P'^^{k) has smaller amplitudes at all scales than 
p2LPTj^^-j is consistent with the result of Crocce, Pueblas & Scoccimarro 2006. This is because 
the higher-order terms in ZA do not necessarily correspond to the growing solutions. The plot 
indicates that differences in the initial conditions affect the power spectrum at the few percent 
level. We confirmed that this difference between 2LPT and ZA reduces to subpercent when 
we increase the starting redshift to = 127. However, we need a very large number of mesh 
points in the calculation of the PM force to avoid large errors introduced by the tree force 
that arise when the density field is nearly uniform. This is very time consuming. Thus we 
need 2LPT initial condition for the purpose of present analysis if we start the simulations at 
Z\^i = 31. We also confirm that our 2LPT result is not affected when we change the starting 
redshift to = 63 or Zi^i = 127. 

A. 1.2. Comparison using different N-body solvers 

Next we compare three different N-body solvers; a Tree-Particle-Mesh (Tree-PM) solver 
Gadget2 (Springel 2005), another Tree-PM code TPM of Bode & Ostriker (2003), and a Particle- 
Particle-Particle-Mesh (P^M) code developed in Jing & Suto (1998, 2002; hereafter JS02). We 
use identical initial conditions as that of one of our fiducial runs, 512^ particles in a box 
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Fig. 8. Wc plot the evolution of tlie power spectrum from two simulations using a ZA initial condition 
and a 2LPT. To see the difference clearly, we show the ratio pZA^p2LPT ^ number of output times. 

of 1000/i~^Mpc. We stop the simulations at z = 3. In addition to the Tree-PM and P^M 
calculations, we also run simulations using only the PM part of the three codes. We measure 
the power spectra for the simulation output at 2 = 3 and correct for the effect of finite volume 
(see section 4 above) so that the comparison is clearer. 

Figure 9 shows the relative difference from the RPT prediction, [-PcOTrccted(^) ~ 
P^^^(A;)]/P'^"(A;). The upper panel shows the results when we include the short-range forces 
(i.e., tree or particle-particle force). The result of TPM, JS02 and Gadget2 are shown by dotted, 
dashed and solid lines with errorbars, respectively. We omit the errorbars for TPM and JS02 
as these are almost the same as for Gadget2, which is given by equation (13). The overall be- 
haviours for the three codes are very consistent, and the difference is important (^ 1%) only at 
smaller scale than k^^^ of RPT/CLA shown by the vertical arrow. In the bottom panel we plot 
the results of the three codes when we use only PM forces. The degree of dispersion among the 
three is larger than that in the upper panel, and is ~ 3% at of RPT/CLA. The short-range 
force is helpful to increase the accuracy at this level at this scale. We conclude that our main 
result is not affected so much by the choice of particular A^-body solver when we turn on the 
short-range force. 
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Fig. 9. A comparison of different iV-body simulation codes. We plot [Pj^j^ertcdl^) " P^^^{k)]/P'''"{k) at 
z = 3. Different line types correspond to TPM (dotted), JS02 (dashed) and Gadget2 (solid), respectively. The 
vertical arrows show the limit wavenmnber fc'™ for RPT/CLA. The errorbars are given by equation (13) for 
Gadget2. Upper: Results of Trce-PM and PP-PM calculations. Lower: Results of PM only calculations. 

A. 1.3. Effects of finite box size and resolution 

The finite box size introduces an additional variance in the power spectrum as shown 
in Takahashi et al. (2008) and this work. It is possible that this finiteness affects not only the 
variance but also the mean value. The resolution of simulations may also affect the results. 
We thus investigate the box size/resolution dependence of the power spectra from TV-body 
simulations. 

In order to increase the spatial resolution, we run simulations varying the box size [1000 
(4 realizations; fiducial), 500 (7 realizations) and 250h~^Mpc (9 realizations)] but keeping the 
number of particles to be the same as our fiducial runs (512^) to check the convergence. For 
this test, we use outputs at z = 3. The results are shown in figure 10. Symbols correspond to 
Lbox = 1000/i~^Mpc (circle), 500/i"^Mpc (triangle) and 250/i~^Mpc (square). We adopt the bin 
width AA; = O.Ol/iMpc"^ in this section for clarity. We also calculate the RPT predictions with 
the effect of finite box size by introducing a cut-off scale (fcbox = Svr/Lbox) in the integration, 
which are shown in the left panel. The solid line is the normal RPT result (without cut-off), 
and the other lines take account of the effect of finite box size (short-dashed: lOOO/i-^Mpc, 
long-dashed: 500/i~^Mpc and dot-dashed: 250/i~^Mpc). Note that we show only the results 
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without cutoff in the text. In the right panel, we plot a comparison with the SPT prediction 
(no cut-off) to show the dependence on the box size of the valid ranges of SPT. 

In the left panel, one may notice the difference of theoretical lines at large wavenumbers. 
This reflects the finiteness of the box: short-dashed line (Lbox = lOOO/i^^ Mpc case) deviates 
from solid line (infinite volume case) only at A; > 0.3h Mpc and is very small (~ 0.5% at 
k = OAh Mpc ~^), whereas long-dashed and dash-dotted ones (Lbox = 500 and 250h~^ Mpc cases) 
diverge at smaller wavenumbers and to a greater degree (~ 1.5% and ~ 3% aX k = OAh Mpc 
respectively). Although the results of A^-body simulations with smaller box sizes have larger 
variances due to smaller total volumes than the fiducial runs, box size dependence appears at 
small scales. The runs with smaller box size tend to have larger amplitudes, which is consistent 
with RPT predictions with cutoff scales. We thus conclude that the results of simulations also 
suffer from a systematic effect due to finite box size, and adopt Lbox = 1000/i~^Mpc for our 
fiducial model, where the effect is ^ 0.5%. 

Next, the convergence regime of SPT seems wider if one believes the results of the smaller 
box size runs in the right panel. We consider these wider valid ranges are just coincidences due 
to the systematic effects in the simulations with smaller box sizes and should not be trusted. 
One should keep this in mind when one uses the results of simulations with small box size. Even 
the Millennium simulation (Springel et al. 2005) has only 500/i~^Mpc on each side, although it 
has many more particles (2160'^) than ours and thus the spacial resolution is better. This box 
size introduces ^ 1% systematics at A; ^ 0.35/iMpc~^ according to our figure 10. Part of the 
disagreement in the valid range of SPT between our results and that of Jeong & Komatsu (2006) 
may be due to this effect, as they use Lbox = 512, 256, 128 and 64/i~^Mpc. 
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Fig. 10. We compare results of TV-body simulations for different box sizes at z = 3. circle: 1000/).~^Mpc, 
triangle: 500/i~^Mpc, square: 250/i~^Mpc. Left: comparison with the RPT prediction (dotted line). We 
also plot predictions of RPT taking account of the finiteness of the volume by introducing a cut off scale: 
solid: short-dashed: lOOO/i-^Mpc (dashed), SOO/i-^Mpc (long-dashcd) and 250/i~iMpc (dot-dashed). 
Right: comparison with the SPT prediction, in which we do not take account of the finiteness effect. 
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